﻿using LxCFD;
using System;
using System.Collections.Generic;

namespace LxCFD
{
    public partial class MatrixBase : MatrixADbBase
    {
        public double[] ALUO;
        public double[] ALUN;
        #region 初始化及公用
        public override void CheckCreate(FvZoneBase zone)
        {
            base.CheckCreate(zone);
            LxTool.CheckCreateArray(ref ALUO, Faces.Length);
            LxTool.CheckCreateArray(ref ALUN, Faces.Length);
        }
        #endregion
        #region 计算
        protected void set_ALU_0(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ALUO[i1] = ALUN[i1] = 0;
            }
        }
        protected double calALUProductValue(ref FvCell cell, double[] fromArray)
        {
            double aluProductValue = 0;
            foreach (int faceIndex in cell.OFaces)
            {
                aluProductValue += ALUO[faceIndex] * fromArray[Faces[faceIndex].NSideCell];
            }
            foreach (int faceIndex in cell.NFaces)
            {
                aluProductValue += ALUN[faceIndex] * fromArray[Faces[faceIndex].OSideCell];
            }
            return aluProductValue;
        }
        protected void cal_x(IndexNumber indexNumber, double[] fromArray, double[] toArray, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                toArray[i1] = (b[i1] - calALUProductValue(ref Cells[i1], fromArray)) / ADiag[i1];
            }
        }
        protected void cal_x_Copy(IndexNumber indexNumber, double[] fromArray, double[] toArray, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                toArray[i1] = fromArray[i1];
            }
        }
        #endregion
    }
    public partial class MatrixBase : ISolvableZone
    {
        protected double Max_Rc;
        protected double Max_acPhic;
        public void CalResidual(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double max_acphic = 0, max_Rc = 0;
            for (; i1 < endp1; i1++)
            {
                double acphic = ADiag[i1] * x[i1];
                double Rc = Math.Abs(acphic + calALUProductValue(ref Cells[i1], x) - b[i1]);
                acphic = Math.Abs(acphic);
                if (acphic > max_acphic) max_acphic = acphic;
                if (Rc > max_Rc) max_Rc = Rc;
            }
            ThreadsMath.KeepMax(ref this.Max_Rc, max_Rc);
            ThreadsMath.KeepMax(ref this.Max_acPhic, max_acphic);
        }
        public void GetResidual_Reset(out double Rc, out double acPhic)
        {
            Rc = this.Max_Rc;
            acPhic = this.Max_acPhic;
            this.Max_Rc = 0;
            this.Max_acPhic = 0;
        }
    }
    public partial class MatrixBase 
    {
        public void CalByJacobi_Step1(int taskIndex)
        {
            foreach (int b in e_FixValue)
            {
                cal_x_Copy(Zone.GetThread(b).ThreadIndex, x, x_Common, taskIndex);
            }
            cal_x(InnerIndex, x, x_Common, taskIndex);
        }
        public void CalByJacobi_Step2(int taskIndex)
        {
            cal_x(InnerIndex, x_Common, x, taskIndex);
        }
    }
    public partial class MatrixBase : IExplicitSolvableZone
    {
        public void CalByExplicit_Step1(int taskIndex)
        {
            cal_x(InnerIndex, x, x_Common, taskIndex);
        }
        public void CalByExplicit_Step2(int taskIndex)
        {
            cal_x_Copy(InnerIndex, x_Common, x, taskIndex);
        }
    }
    public partial class MatrixBase : IAMGSolvableZone
    {
        public void Restrict_Ar(out double A, ref MGCell mgCell)
        {
            double Atemp = 0;
            double rtemp = 0;
            foreach (int finecell in mgCell.FineCells)
            {
                double atemp = ADiag[finecell];
                Atemp += atemp;
                rtemp += b[finecell] - calALUProductValue(ref Cells[finecell], x) - atemp * x[finecell];
            }
            foreach (int innerfineface in mgCell.InnerFineFaces)
            {
                Atemp += ALUO[innerfineface] + ALUN[innerfineface];
            }
            A = Atemp;
            mgCell.r = rtemp;
        }
        public void Restrict_ALU(out double A_lu_o, out double A_lu_n, ref MGFace coarseFace)
        {
            double O2N = 0;
            double N2O = 0;
            foreach (int fineface in coarseFace.FineFaces_O2N)
            {
                O2N += ALUO[fineface];
                N2O += ALUN[fineface];
            }
            foreach (int fineface in coarseFace.FineFaces_N2O)
            {
                O2N += ALUN[fineface];
                N2O += ALUO[fineface];
            }
            A_lu_o = O2N;
            A_lu_n = N2O;
        }
        public void Restrict_r(ref MGCell coarseCell)
        {
            double rtemp = 0;
            foreach (int finecell in coarseCell.FineCells)
            {
                rtemp += b[finecell] - calALUProductValue(ref Cells[finecell], x) - ADiag[finecell] * x[finecell];
            }
            coarseCell.r = rtemp;
        }
        public double[] GetProlongArray()
        {
            return x;
        }
        public void GetBounds(out List<int> e_fixValue, out List<int> e_free, out List<int> e_couple)
        {
            e_fixValue = this.e_FixValue;
            e_free = this.e_ALUO0;
            e_couple = this.e_Couple;
        }
    }
    public partial class MatrixBase
    {
        public void Initialize_Step1_r0rp1rho(int taskIndex)
        {
            var r0 = Zone.BiCGSolver.r0;
            var r = Zone.BiCGSolver.r;
            var p = Zone.BiCGSolver.p;
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double temp = 0;
            for (; i1 < endp1; i1++)
            {
                double rt = r[i1] = r0[i1] = p[i1] = b[i1] - calALUProductValue(ref Cells[i1], x) - ADiag[i1] * x[i1];
                temp += rt * rt;
            }
            ThreadsMath.AddToSum(temp);
            foreach (int b in e_FixValue)
            {
                Zone.BiCGSolver.Initialize0ValueInArray(Zone.GetThread(b).ThreadIndex, taskIndex);
            }
        }
        public void Cal_rhoPart1(int taskIndex)
        {
            Zone.BiCGSolver.Cal_r_Dot_r0(InnerIndex, taskIndex);
        }
        public void Cal_pArray(int taskIndex)
        {
            Zone.BiCGSolver.Cal_pArray(InnerIndex, taskIndex);
        }
        public void Cal_yArray(int taskIndex)
        {
            Zone.BiCGSolver.Cal_yArray(InnerIndex, ADiag, taskIndex);
        }
        public void Cal_vArray_alpha1(int taskIndex)
        {
            double[] v = Zone.BiCGSolver.v;
            double[] y = Zone.BiCGSolver.y;
            double[] r0 = Zone.BiCGSolver.r0;
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double temp = 0;
            for (; i1 < endp1; i1++)
            {
                double vt = v[i1] = ADiag[i1] * y[i1] + calALUProductValue(ref Cells[i1], y);
                temp += r0[i1] * vt;
            }
            ThreadsMath.AddToSum(temp);
        }
        public void Cal_sArray(int taskIndex)
        {
            Zone.BiCGSolver.Cal_sArray(InnerIndex, taskIndex);
        }
        public void Cal_zArray(int taskIndex)
        {
            Zone.BiCGSolver.Cal_zArray(InnerIndex, ADiag, taskIndex);
        }
        public void Cal_tArray_omega1(int taskIndex)
        {
            double[] s = Zone.BiCGSolver.r;
            double[] z = Zone.BiCGSolver.z;
            double[] t = Zone.BiCGSolver.t;
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double temp1 = 0, temp2 = 0;
            for (; i1 < endp1; i1++)
            {
                double st = s[i1];
                double tt = t[i1] = ADiag[i1] * z[i1] + calALUProductValue(ref Cells[i1], z);
                temp1 += tt * st;
                temp2 += tt * tt;
            }
            ThreadsMath.AddToSum(temp1, temp2);
        }
        public void Cal_xArray(int taskIndex)
        {
            Zone.BiCGSolver.Cal_xArray(InnerIndex, x, taskIndex);
        }
        public void Cal_rArray(int taskIndex)
        {
            Zone.BiCGSolver.Cal_rArray(InnerIndex, taskIndex);
        }
    }
}
//public void SendRestrict_A_LU_O(out double A_lu_o, ref MGFace coarseFace)
//{
//    double O2N = 0;
//    foreach (int fineface in coarseFace.FineFaces_O2N)
//    {
//        O2N += ALUO[fineface];
//    }
//    foreach (int fineface in coarseFace.FineFaces_N2O)
//    {
//        O2N += ALUN[fineface];
//    }
//    A_lu_o = O2N;
//}
//public void CalByGS_Forward(int taskIndex)
//{
//    cal_x(InnerIndex, x, x, taskIndex);
//}
//public void CalByGS_Backward(int taskIndex)
//{
//    InnerIndex.CellIDs.GetSeg(out int start, out int i2, taskIndex);
//    i2--;
//    for (; i2 >= start; i2--)
//    {
//        x[i2] = (b[i2] - calALUProductValue(ref cells[i2], x)) / ADiag[i2];
//    }
//}
//if (AppRes.computeCenter.isTransient)
//{
//    if (transientRecorder == null)
//    {
//        transientRecorder = new TransientDataRecorder();
//        transientRecorder.Create(this);
//    }

//}
//else
//{
//    transientRecorder = null;
//}